librarian::shelf(clusterProfiler,
                 enrichplot,
                 tidyverse,
                 ggpubr,
                 patchwork)
source('00_util_scripts/mod_bplot.R')

enri_simplify <- clusterProfiler::simplify

himrsa_fc <- read_csv('mission/FPP/zww_sa_mice/results/himrsa-pbs-logFC-edgeR-QLF.csv')

up24 <- himrsa_fc |>
  filter(str_detect(test, '24h') & FDR < .05 & logFC > 2)

down24 <- himrsa_fc |>
  filter(str_detect(test, '24h') & FDR < .05 & logFC < -2)

up8 <- himrsa_fc |>
  filter(str_detect(test, '8h') & FDR < .05 & logFC > 2)

down8 <- himrsa_fc |>
  filter(str_detect(test, '8h') & FDR < .05 & logFC < -2)

# 8h ----------
### FIG: 8h upreg gobp =============
up8_ora_go <- up8 |>
  pull(gene_symbol) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = "BP",
           readable = TRUE)

up8_ora_go@result |>
  head(10) |>
  publish_enrichment()

publish_pdf('micebulk/micebulk_8h_upgo.pdf', width = 65)

## downreg ----------
down8_ora_go <- down8 |>
  pull(gene_symbol) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = "BP",
           readable = TRUE)

down8_ora_go@result |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway 8h after infection')

ggsave('mice_bulk_8h_downgo.pdf',
       width = 7, height = 4)

## GSEA ---------------
gse.8h <- himrsa_fc |>
  filter(str_detect(test, '8h') & FDR < .05) |>
  pull(logFC, name = gene_symbol) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', 'org.Mm.eg.db', 'SYMBOL')

gse.8h@result |>
  write_csv('mission/FPP/zww_sa_mice/results/bulk.8h.gsea.csv')

gse.8h.tb <- read_csv('mission/FPP/zww_sa_mice/results/bulk.8h.gsea.csv')

gse.8h.tb |>
  filter(str_detect(Description, 'T')) |>
  DT::datatable()


# 24h ----------
### FIG: 24h upreg gobp =============
up24_ora_go <- up24 |>
  pull(gene_symbol) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = "BP",
           readable = TRUE) |>
  enri_simplify()

up24_ora_go@result |>
  mutate(p.adjust = signif(p.adjust, 2)) |>
  slice_max(Count, n=10) |>
  publish_enrichment()

publish_pdf('micebulk/micebulk_24h_upgo.pdf', width = 65)

## downreg
down24_ora_go <- down24 |>
  pull(gene_symbol) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = "BP",
           readable = TRUE) |>
  enri_simplify()

down24_ora_go@result |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway 24h after infection')

ggsave('mice_bulk_24h_downgo.pdf',
       width = 7, height = 4)

## GSEA ---------------
gse.24h <- himrsa_fc |>
  filter(str_detect(test, '24h') & FDR < .05) |>
  pull(logFC, name = gene_symbol) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', 'org.Mm.eg.db', 'SYMBOL') |>
  simpl_enrich()

gse.24h@result |>
  write_csv('mission/FPP/zww_sa_mice/results/bulk.24h.gsea.csv')

gse.24h.tb <- read_csv('mission/FPP/zww_sa_mice/results/bulk.24h.gsea.csv')

gse.24h.tb |>
  filter(str_detect(Description, 'T')) |>
  DT::datatable()

# MVK kegg --------
gene_of_mva <- c('Hmgcr','Hmgcs1','Hmgcs2','Fdps','Mvd','Idi1','Idi2')

kegg_mva <- c(gene_of_mva, 'Pmvk','Acat1','Acat2','Mvk')

himrsa_fc |>
  filter(gene_symbol %in% kegg_mva) |>
  mutate(test = fct_relevel(test, '8h vs PBS'),
         type = case_when(FDR < .05 & logFC > 0 ~ 'Upregulated',
                          FDR < .05 & logFC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(gene_symbol, logFC, fill = type)) +
  geom_col() +
  facet_wrap(~test) +
  scale_fill_manual(values = c('blue','grey','red')) +
  coord_flip() +
  theme_pubr()

ggsave('mice_bulk_mva_logfc_barplot.pdf',
       width = 6, height = 4)
